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This handout provides the essential elements needed to understand finite random matrix theory. A 
cursory observation should reveal that the tools for infinite random matrix theory are quite different from 
the tools for finite random matrix theory. Nonetheless, there are significantly more published applications 
that use finite random matrix theory as opposed to infinite random matrix theory. Our belief is that many 
of the results that have been historically derived using finite random matrix theory can be reformulated and 
answered using infinite random matrix theory. In this sense, it is worth recognizing that in many applications 
it is an integral of a function of the eigenvalues that is more important that the mere distribution of the 
eigenvalues. For finite random matrix theory, the tools that often come into play when setting up such 
integrals are the Matrix Jacobians, the Joint Eigenvalue Densities and the Cauchy-Binet theorem. We 
describe these in subsequent sections. 



1 Matrix and Vector Differentiation 

In this section, we concern ourselves with the differentiation of matrices. Differentiating matrix and vector 
functions is not significantly harder than differentiating scalar functions except that we need notation to 
keep track of the variables. We titled this section "matrix and vector" differentiation, but of course it is 
the function that we differentiate. The matrix or vector is just a notational package for the scalar functions 
involved. In the end, a derivative is nothing more than the "linearization" of all the involved functions. 
We find it useful to think of this linearization both symbolically (for manipulative purposes) as well as 
numerically (in the sense of small numerical perturbations). The differential notation idea captures these 
viewpoints very well. 

We begin with the familiar product rule for scalars, 

d{uv) ~ M(di)) + v{du), 

from which we can derive that d(a;'^) = 3x^dx. We refer to da: as a differential. 

We all unconsciously interpret the "dx" symbolically as well as numerically. Sometimes it is nice to 
confirm on a computer that 

^^i±4— ^«3x^ (1) 

I do this by taking x to be 1 or 2 or rsindnd) and e to be .001 or .0001. 
The product rule holds for matrices as well: 

d{UV) = U{dV) + {dU)V. 

In the examples we will see some symbolic and numerical interpretations. 

Example 1: Y ^ 

We use the product rule to differentiate Y{X) = to obtain that 



d{X^) = X^{dX) + X{dX)X + {dX)X^ . 



When I introduce undergraduate students to matrix multiplication, I tell them that matrices are like scalars, 
except that they do not commute. 

The numerical (or first order perturbation theory) interpretation applies, but it may seem less familiar 
at first. Numerically take X=randn(n) and E=randn(n) for e = .001 say, and then compute 

This is the matrix version of (1). Holding X fixed and allowing E to vary, the right-hand side is a linear 
function of E. There is no simpler form possible. 

Symbolically (or numerically) one can take AX = Eki which is the matrix that has a one in element (/c, I) 
and elsewhere. Then wc can write down the matrix of partial derivatives: 

gv3 

—— = X^iEki) + X{Eki)X + {Eki)X^ . 
oxki 

As we let fc, I vary over all possible indices, we obtain all the information we need to compute the derivative 
in any general direction E. 

In general, the directional derivative of Yij{X) in the direction dX is given by {dY)ij. For a particular 
matrix X, dY{X) is a matrix of directional derivatives corresponding to a first order perturbation in the 
direction E ~ dX. It is a matrix of linear functions corresponding to the linearization of Y{X) about X. 

Structured Perturbations 

Wc sometimes restrict our _B to be a structured perturbation. For example if X is triangular, symmetric, 
antisymmetric, or even sparse then often we wish to restrict E so that the pattern is maintained in the 
perturbed matrix as well. An important case occurs when X is orthogonal. We will see in an example below 
that wc will want to restrict E so that X'^E is antisymmetric when X is orthogonal. 

Example 2: y = x'^x 

Here y is a scalar and dot products commute so that dy ~ 2x'^dx. When y = 1, x is on the unit sphere. 
To stay on the sphere, we need dy = so that x'^dx = 0, i.e., the tangent to the sphere is perpendicular to 
the sphere. Note the two uses of dy. In the first case it is the change to the squared length of x. In the 
second case, by setting dy = 0, wc find perturbations to x which to first order do not change the length at 
all. Indeed if one computes {x + dx)'^ {x + dx) for a small finite dx, one sees that if x'^dx = then the length 
changes only to second order. Geometrically, one can draw a tangent to a circle. The distance to the circle 
is second order in the distance along the tangent. 

Example 3: y = x'^Ax 

Again y is scalar. We have dy = dx'^Ax + x^Adx. If A is symmetric then dy = 2x'^Adx. 

Example 4: Y = X-^ 

We have that XY ^ I so that X{dY) + idX)Y = so that dY = -X-^dXX-\ 

We recommend that the reader compute e^^{{X + eE)^^ — X^^) numerically and verify that it is equal 
to -X-^EX^'^. 
In other words, 

{X + eE)-^ = X-'^ - eX-^EX-^ + 0{e^) . 

Example 5: / = Q'^Q 

If Q is orthogonal we have that Q^dQ + dQ^Q = so that Q^dQ is antisymmetric. 
In general, d{Q'^Q) ~ Q'^dQ+dQ^^Q, but with no orthogonality condition on Q, there is no anti-symmetry 
condition on Q^dQ. 

If y is a scalar function oi xi, X2, ■ ■ ■ , Xn then we have the "chain rule" 



Often we wish to apply the chain rule to every element of a vector or matrix. 

„2 



Let X 
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and Y = X^ 



p + qr 
pr + rs 



pq^ 
qs 



rs 

„2 



then 



dY = XdX + dXX . 



(3) 



2 Matrix Jacobians (getting started) 
2.1 Definition 

Let X G M" and y = y{x) E R" be a difFerentiable fmrction of It is well known that the Jacobian matrix 



J = 



dxi dXn 

dyn dy„ 
\ dxi dXn ' 



dxi 



ij = l,2,. 



evaluated at a point x approximates y{x) by a linear function. Intuitively y{x + 5x) w y{x) + J5x^ i.e., J 
is the matrix that allows us to invoke perturbation theory. The function y may be viewed as performing a 
change of variables. 

Furthermore (intuitively) if a little box of n dimensional volume e surrounds x, then it is transformed by 
y into a parallelopipcd of volume | det J\e around y{x). Therefore the Jacobian \ det J| is the magnification 
factor for volumes. 

If we are integrating some function of y G M" as in J p{y)dy, (where dy = dyi . . - dyn), then when we 
change variables from y to x where y = y{x), then the integral becomes j p(y(x)) |det(^^)| da;. For many 
people this becomes a matter of notation, but one should understand intuitively that the j'acobian tells you 
how volume elements scale. 

The determinant is exactly where the change of variables breaks down. It is always instructive to see 
when this happens. Either there is no "x" locally for each "y" or there are many as in the example of polar 
coordinates at the origin. 

Notation: throughout this book, J denotes the Jacobian matrix. (Sometimes called the derivative or 
simply the Jacobian in the literature.) We will consistently write det J for the Jacobian determinant (un- 
fortunately also called the Jacobian in the literature.) When we say Jacobian, we will be talking about 
both. 

2.2 Simple Examples (n=2) 

We get our feet wet with some simple 2x2 examples. Every reader is familiar with changing scalar variables 
as in 

..2\ 



fix)dx^ J f{y'){2y)dy. 

We want the reader to be just as comfortable when / is a scalar function of a matrix and we change X = Y'^: 

f{X){dX)^ I /(y 2) ( Jacobian) (dr). 



One can compute all of the 2 by 2 Jacobians that follow by hand, but in some cases it can be tedious 
and hard to get right on the first try. Code 8.1 in MATLAB takes away the drudgery and gives the right 
answer. Later we will learn fancy ways to get the answer without too much drudgery and also without the 
aid of a computer. We now consider the following examples: 



2x2 Example 1: Matrix Square {Y = X'^) 



2x2 Example 2 
2x2 Example 3 
2x2 Example 4 
2x2 Example 5 
2x2 Example 6 
2x2 Example 7 
2x2 Example 8 
2x2 Example 9 



Matrix Cube {Y = X^) 

Matrix Inverse {Y = X-'^) 

Linear Transformation (Y = AX + B) 

The LU Decomposition {X = LU) 

A Symmetric Decomposition (5* — DMD) 

Traceless Symmetric — Polar Decomposition {S = QAQ^, trS* = 0) 
The Symmetric Eigenvalue Problem {S — QhQ^) 
Symmetric Congruence {Y = A^SA) 



Discussion: 

Example 1: Matrix Square (F = X^) 
P q 



With X 



and Y = X the Jacobian matrix of interest is 
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J = 



On this first example we label the columns and rows so that the elements correspond to the definition 
J = (^^^y Later we will omit the labels. We invite readers to compare with Equation (3). We see that 
the Jacobian matrix and the differential contain the same information. We can compute then 

det J = 4(p + s)^{sp - qr) = 4(trX)2 det(X) . 

Notice that breakdown occurs if X is singular or has trace zero. 

Example 2: Matrix Cube (Y = X^) 
p q 



With X = 



and y = 



J = 



2rp + sr 

„2 



3p- + 2qr pq + q{p + s) 

2rp + sr p'^ + 2 qr + (p + s) s r" 

2pq + qs q^ p{p + s) + 2qr + 

qr pq + 2qs r (p + s) + sr 



qr 
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2 gr + 3 



so that 



det J = 9(sp - qrf{qr+p^ + s^ + spf ^ -(det X)2(tr + {tr X 



i2\2 



Breakdown occurs if X is singular or if the eigenvalue ratio is a complex cube root of unity. 

Example 3: Matrix Inverse {Y = X^^) 
p q 



With X = 



r s 



and Y = X- 
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so that 

Breakdown occurs if X is singular. 



det J= (detX)-'' 



Example 4: Linear Transformation {Y = AX + B) 
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The Jacobian matrix has two copies of the constant matrix A so that det J = det A^ — (det A)^. Breakdown 
occurs if A is singular. 

Example 5: The LU Decomposition {X = LU) 

The LU Decomposition computes a lower triangular L with ones on the diagonal and an upper triangular 
U such that X = LU. 

For a general 2x2 matrix it takes the form 
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which exists when 0. 

Notice the function of four variables might be written: 



2/1 
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2/3 
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The Jacobian matrix is itself lower triangular 
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so that det J = I /p. Breakdown occurs when p ~ 0. 



Example 6: A Symmetric Decomposition (S = DMD) 



Any symmetric matrix X ~ 



p r 
r s 



may be written as 



X = DMD where D 



VP 





and M 



1 

rjJps 1 



Since X is symmetric, there arc three independent variables, and thus the Jacobian matrix is 3 x 3. The 
three independent elements in D and M may be thought of as functions of p, r, and s: namely 

2/1 = Vp 

y2 = \fs and 

2/3 = rj ^ps 



The Jacobian matrix is 



7f « 



so that 



Breakdown occurs if p or s is 0. 



detJ 



r/s 2 
1 



Example 7: Traceless Symmetric = Polar Decomposition {S = QAQ'^, trS* = 0) 

The reader will recall the usual definition of polar coordinates. If (p, s) are Cartesian coordinates, then 
the angle is 6 ~ arctan(s/p) and the radius is r = ^yp^+s^. If we take a symmetric traceless 2x2 matrix 

S 



p s 
s —p 



and compute its eigenvalue and eigenvector decomposition, we find that the eigendecomposition is mathe- 
matically equivalent to the familiar transformation between Cartesian and polar coordinates. Indeed one of 
the eigenvectors of S is (cos 0, sin 0), where cj) = 9/2. The Jacobian matrix is 



J = 



-S P 



The Jacobian is the inverse of the radius. This corresponds to the familiar formula using the more usual 
notation dxdy/r — drd9 so that dct J ~ l/r. Breakdown occurs when r ~ 0. 

Example 8: The Symmetric Eigenvalue Problem {S = QKQ^) 

We compute the Jacobian for the general symmetric eigenvalue and eigenvector decomposition. Let 
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We can compute the eigenvectors and eigenvalues of S directly in MATLAB and compute the Jacobian of 
the two eigenvalues and the eigenvector angle, but when we tried with the Maple toolbox we found that the 
symbolic toolbox did not handle "arctan" very well. Instead we found it easy to compute the Jacobian in 
the other direction. 

We write S = QAQ^ , where Q is 2 x 2 orthogonal and A is diagonal. The Jacobian is 
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so that 

det J = Ai 

Breakdown occurs if 5 is a multiple of the identity. 



Example 9: Symmetric Congruence {Y = A^SA) 



Let Y = A^SA, where Y and S are symmetric, but A ~ 



J = 



a 

fe2 



c 



2 ca 
2db 



ab cd cb + ad 



is arbitrary. The Jacobian matrix is 



and det J = (det^)^. 

The cube on the determinant tends to surprise many people. Can you guess what it is for an n x n 
symmetric matrix (1^ = A'^SA)? The answer (det J = (det yl)""'"^)is in Example 3 of Section 3. 



%jacobiaii2by2 .m 

"/oCode 8.1 of Random Eigenvalues by Alan Edelman 

"/oExperiment : Compute the Jacobian of a 2x2 matrix function 

"/Comment: Symbolic tools are not perfect. The author 

7„ exercised care in choosing the variables. 

syms pqrsabcdtele2 
X=[p q ; r s] ; A= [a b;c d] ; 

%7o Compute Jacobians 



)), JAC_square =factor(det(J)) 
)), JAC_cube =factor(det(J)) 



Y=X~2; J= jacobian (Y(:),X(: 

Y=X~3; J= jacobian (Y(:),X(: 

Y=inv(X); J=jacobian(Y( : ) ,X( : 

Y=A*X; J= jacobian (Y(:),X(: 

Y=[p q;r/p det(X)/p]; J=jacobian(Y( : ) ,X( : ) ) , JAC_lu =f actor (det (J) ) 



)), JAC_inv =factor(det(J)) 
)), JAC_linear =factor(det(J)) 



x=[p s r] ;y=[sqrt(p) sqrt(s) r/(sqrt(p)*sqrt(s))] ; 

J=jacobian(y,x) , JAC_DMD =factor(det(J)) 

x= [p s] ; y= [ sqrt (p~2+s~2) atan(s/p)] ; 

J=jacobian(y ,x) , JAC_notrace =factor(det(J)) 

Q=[cos(t) -sin(t) ; sin(t) cos(t)]; 
D=[el 0;0 e2] ; Y=q*D*Q . ' ; 

y=[Y(l,l) Y(2,2) Yd, 2)]; x= [t el e2] ; 

J=jacobian(y ,x) , JAC_symeig =simplif y (det ( J) ) 
X=[p s;s r] ; Y=A.'*X*A; 
y=[Y(l,l) Y(2,2) Yd, 2)]; x= [p r s] ; 

J=jacobian(y ,x) , JAC_symcong =factor(det(J)) 



Code 1 



3 Y = BXA^ and the Kronecker Product 



3.1 Jacobian of y = BXA^ (Kronecker Product Approach) 

There is a "nuts and bolts" approach to calculate some Jacobian determinants. A good example function is 
the matrix inverse Y — X^^. We recall from Example 4 of Section 1 that 

dY = ~X-^dXX~^ . 

In words, the perturbation dX is multiplied on the left and on the right by a fixed matrix. When this 
happens we are in a "Kronecker Product" situation, and can instantly write down the Jacobian. 

We provide two definitions of the Kronecker Product for square matrices A e M"-" to B e E'"-"' . 

See [1] for a nice discussion of Kronecker products. 

Operator Definition 

A(g)B is the operator from X e M™'" to y e M'"'" where Y = BXA'^. We write 



{A®B)X = BXA^ . 

Matrix Definition (Tensor Product) 
A(E) B is the matrix 

^ aiiB . . . aim^B \ 
A®B= ■. : . (4) 

\ ^milB . . . Ciinim2B j 

The following theorem is important for applications. 
Theorem 1. d(±{A®B) = (det yl)™(det B)" 

Application: If y = then dY = -(A"^ ® A"!) dA so that 

|det J| = |detA|-2". 

Notational Note: The correspondence between the operator definition and matrix definition is worth 
spelling out. It corresponds to the following identity in MATLAB 

Y = B * X * A' 

Y(:) = kron(A,B) * X(:) 

"/o The second line does not change Y 

Here kron(A,B) is exactly the matrix in Equation (4), and X(:) is the column vector consisting of the 
columns of A stacked on top of each other. (In computer science this is known as storing an array in 
"column major" order.) Many authors write vec(X), where we use X(:). Concretely, we have that 

vec {BXA^ = (A ® B)vcc(A) 

where A (g) B is as in (4) . 
Proofs of Theorem 8.1: 

Assume A and B are diagonalizable, with Aui = XiUi (i = 1, . . . ,ri) and Bvi ~ i^iiVi {i = 1, . . . ,m). 
Let Alij — ViuJ . The mn matrices My form a basis for R™" and they are eigenmatrices of our map since 
BAIijA^ = jiiXjAIij. The determinant is 

W ti.,Xj or (det^)"(detS)". (5) 

l<i<Ti 
l<j<m 



The assumption of diagonalizability is not important. 



Also one can directly manipulate the matrix since obviously 

A®B^ {A(g)I){I®B) 

as operators, and det I ® B = (det i?)" and det A® I which can be reshuffled into det I ® A = (det A)"^ . 

Other proofs may be obtained by working with the "LC/" decomposition of A and B. the SVD of A and 
B, the QR decomposition, and many others. 

Mathematical Note: That the operator Y = BXA^ can be expressed as a matrix is a consequence of 
linearity: 

B{ciXi + C2X2)A^ = ciBXiA^ + C2BX2A^ 
i.e. {A®B){ciXi + c2X2) = ci{A® B)Xi + C2{A® B)X2 

It is important to realize that a linear transformation from R™'" to M™-" is defined by an element of 
M™"'™". i.e., by the m^n^ entries of an mn x mn matrix. The transformation defined by Kronecker products 
is an + subspace of this m^n^ dimensional space. 

Some Kronecker product properties: 

1. {A®B)^=A'^(g>B^ 

2. {A(^B)-^ = 

3. det(yl ® B) ^ det(A)™ det(B)", A e M"'", B e M™''" 

4. ti{A(g)B) = tr(A)tr(B) 

5. A (g) i? is orthogonal if A and S are orthogonal 

6. {A (g) B)(C ® Z?) = (AC) ® {BD) 

7. If Au = Am, and i?w = ^^v, then if X = wu^, then BXA'^= XfxX, and also AX'^B^ X^X'^ . Therefore 
A® B and B A have the same eigenvalues, and transposed eigenvectors. 

It is easy to see that property 5 holds, since if A and B are orthogonal {A (g B)'^ = A'^® B^ = A~^ (g) B~^ = 

{A(g}B)-\ 

Linear Subspace Kronecker Products 

Some researchers consider the symmetric Kronecker product (ggym- In fact it has become clear that one 
might consider an anti-symmetric, upper-triangular or even a Toeplitz Kronecker product. We formulate a 
general notion: 

Definition: Let S denote a linear subspace of R™" and tt^ a projection onto S. If A e R"-" and B E 
then we define {A g) B)sX ^ tts{BXA^ for X G 5. 

Comments 

1. If 5 is the set of symmetric matrices then m — n and 

BXA'^+AXB^ 
(A (E) B)^y^X = 



2. If iS is the set of anti-sym matrices, then m — n 



BXA^+AXB^ 
{A g) Bjanti-'^ = 7, as well 



but this matrix is anti-sym. 



3. If S is the set of upper triangular matrices, then m = n and [A (g) i?)uppor = upper(_BX^"^. 

4. We recall that is a projection if tt^ is linear and for all X E K"'", tt^X e S. 

5. Jacobian of ( 



vsym ■ 



{A ® B%yn,X 



BXA^+AXB'^ 



If B = A then 



det J = J^AiAj = (det^)"+^ 



by considering eigenmatrices E = UiuJ for i < j. 

6. Jacobian of iKiuppcr: 

Special case: A lower triangular, B upper triangular so that 



(A ® B)„.X = BXA^ 



since BXA^ is upper triangular. 



The eigenvalues of A and B are Xi = An and = Bjj respectively, where Aui = XiUi and Bui = fiiVi 
{i = 1, . . . ,n). The matrix Alij = ViuJ for « < j is upper triangular since Vi and Uj are zero below the 
ith and above the jth component respectively, ("^he eigenvectors of a triangular matrix are triangular.) 



M,,. 




for i < j . 



Since the Mij are a basis for upper triangular matrices BAIijA'^ — jiiXjMij. We then have 

detJ = nM.A, = {\,\l\l...\l){^i-,^4~^^ir^...^lr.) 

i<j 

= (^11^22^33 . . . An7i){BllB22^ B^^j^'^ . . . B„„) . 

Note that J is singular if an only if ^ or i? is. 

7. Jacobian of (8)Toopiitz 

Let X be a Toeplitz matrix. We can define 

{A (^Toeplitz B)X = Toeplitz(SX^^ 

where Toeplitz averages every diagonal. 



4 Jacobians of Linear Functions, Powers and Inverses 

The Jacobian of a linear map is just the determinant. This determinant is not always easily computed. The 
dimension of the underlying space of matrices plays a role. For example the Jacobian of F = 2X is 2" for 
X G M"'^", 2^2^' for upper triangular or symmetric X, 2~t ' for antisymmetric X, and 2" for diagonal 
X. 

We will concentrate on general real matrices X and explore the symmetric case and triangular case as 
well when appropriate. 



4.0.1 General Matrices 

Let Y ^X^ then 

dY = XdX + dXX or dY = I®X + X'^®I 

Using the Eij as before Equation (5), we see that the eigenvalues oi I ® X + X'^ ® I are Ai + \j so that 
det J = ni<,j<„(-^» + Aj) = 2"detxn,<j(A^ + Xjf. When n = 2, we obtain 4detX(trX)^ as seen in 
example 8.2.2. 

For general powers, Y = X^, {p = 1, 2, 3, . . .) 

dY = ^XMXXP-i-*^ 

k=Q 

and det J = [| ("^ Af Af '"^^ = p"(dct X)^-! [] 

l<ij<n \k=0 / i<j 

Of course by computing the Jacobian determinant explicitly in terms of the matrix elements, we can obtain 
an equivalent expression in terms of the elements of X rather than the A^. 

Formally, if wc plug in p = —1. wc find that we can compute det J for Y = X~^ yielding the answer 
given after Theorem 8.1. Indeed one can take any scalar function y = f{x) such as x^, e^, log a;, sin a;, etc. 
The correct formula is 

det J = dct(/'(X)) -Jl 

i<j 

We can prove these using the tools provided by wedge products. 

Often one thinks that to define Y = f{X) for matrices, one needs a power series, but this is not so. If 
X = ZAZ-i, then we can define f{X) = Zf{A)Z-^, where /(A) = diag(/(Ai), . . . ,/(A„)). The formula 
above is correct at any matrix X such that / is differentiable at the eigenvalues. 

When X is of low rank, we can consider the Moore-Penrose Inverse of X. The Jacobian was calculated 
in [2] and [3]. There is a nonsymmctric and symmetric version. 

Exercise. Compute the Jacobian determinant for 

Y = X^/^ and Y^X-^^^. 

Advanced exercise: compute the Jacobian matrix for Y = X^l"^ as an inverse of a sum of two Kronecker 
products. Numerically, a "Sylvester equation" solver is needed. To obtain the Jacobian matrix for Y = X~^/'^ 
one can use the chain rule, given the answer to the previous problem. It is interesting that this is as hard as 
it is. 

4.0.2 Symmetric Matrices 

If X is symmetric, we have that det J = (det Ar)^("+^'' for the map Y = X~^. 

5 Jacobians of Matrix Factorizations (without wedge products) 

Let A be an n X n matrix. In elementary linear algebra, we learn about Gaussian elimination, Gram- 
Schmidt orthogonalization, and the eigendecomposition. All of these ideas may be written compactly as 
matrix factorizations, which in turn may be thought of as a change of variables: 



A, — A, 



f{\)-f{X,) 
A,: — A," 



n/'(A^)n 



Aj - A, 







parameter count 


Gaussian Elimination: 


A= L ■ U 

T T 

unit lower upper 
triangular triangular 


7i(n- \)/2 + n{n + l)/2 


Gram-Schmidt : 


A= Q ■ R 

T T 

Orthogonal upper 

triangular 


n{n- l)/2 + n(n+ l)/2 


Eigenvalue Decomposition: 


A= X ■ A • X-^ 

T T 

eigenvectors eigenvalues 


[n^ — ?i) + 71 

T T 

eigenvector eigenvalue 



Each of these factorizations is a change of variables. Somehow the parameters in A are transformed into 
v? other parameters, though it may not be immediately obvious what these parameters are (the n{n — l)/2 
parameters in Q for example). 

Our goal is to derive the Jacobians for those matrix factorizations. 

5.1 Jacobian of Gaussian Elimination {A = LU) 

In numerical linear algebra texts it is often observed that the memory used to store ^ on a computer may 
be overwritten with the n? variables in L and U . Graphically the n x n matrix A 







^-i_ ^ 


A 













Indeed the same is true for other factorizations. This is not just a happy coincidence, but deeply related to 
the fact that parameters ought not need more storage. 

Theorem 2. If A ~ LU , the Jacobian of the change of variables is 

n 

det J = <r^U22"^ • ■ • u„_i,„_i = 

i=l 



Proof 1: Let A = LU , then using 

dA = LdU + dLU 

= L{dUU-^ + L^^dL)U 

^ {U^® L) ((C/^®uppcr ly^dU + (/ Slower L^^ dL) 

The mapping dU dUU^^ only affects the upper triangular part. Similarly dL L^^dL for the lower 
triangular. We might think of the map in block format: 



dA = U'^®L 



(C/^®uppcr/) \^ ( '^^ 

I Slower L J V dL 



Since the Jacobian determinants of U'^ ® L \sY[ "f^u a-^d of {U'^(E) I) ^ is H ""i/ (and that oi I (E) L is 1), we 
have that det J J] "ir'- ^ 

That is the fancy slick proof. Here is the direct proof. It is of value to see both. 



Proof 2: Let M;, = I * ^ '^^ as in the diaffram above. Since A = LU, we have that 

i-l 

Aij = ^ M^kMkj + for i < j 



and 



Therefore 



fc=i 

^ij = ^ MikMkj + kjUjj for i > j . 
fc=i 

_ r 1 if i<j 



dMio \ Ujj if i > j 



(6) 



Notice that Aij never depends on Mpq when p > i or q > j. Therefore if we order the variables first by row 
and next by column, we see that the Jacobian matrix is lower triangular with diagonal entries given in (6). 
Remembering that the determinant of a lower triangular matrix is the product of the diagonal entries, the 
theorem follows. □ 

Perhaps it should not have been surprising that the condition that A admits a unique LU decomposition 
may be thought of in terms of the Jacobian. Given a matrix A, let pk denote the determinant of the upper 
left fc by A; principal minor. It may be readily verified that un . . . Ukk = Pk and hence the Jacobian is 
P1P2 ■ ■ -Pn-i- The condition that A admits a unique LU decomposition is well known to be that all the 
upper-left principal minors of dimension smaller than n are non-singular. Otherwise the matrix may have 
no LU decomposition. For example. 



1 

1 1 



has no LU factorization. 



It may also have many LU decompositions as does the zero matrix when n > 1. This degeneracy explains the 
need for pivoting strategies (i.e., strategies that reorder the matrix rows and/or columns) for solving linear 
systems of equations even if the computation is done in exact arithmetic. Modern Gaussian elimination 
software for solving linear systems of equations in finite precision include pivoting strategies designed to 
avoid being near a matrix with such a degeneracy. 

Exercise. If A is symmetric positive definite, it may be factored A = LL^ . This is the famous Cholesky 
decomposition of A. How many independent variables are in A? In L? Prove that the Jacobian of the change 
of variables is 



detJ = 2"n^r'" 



i=l 



Research Question. It seems that the existence of a finite algorithm for a matrix factorization is linked 
to a triangular Jacobian matrix. After all, the latter implies only one new variable need be substituted at a 
time. This is the essential idea of Doolittle's or Crout's matrix factorization schemes for LU . 



5.2 Jacobian of Gram-Schmidt (A = QR) 

It is well known that any nonsingular square matrix A may be factored uniquely into orthogonal times upper 
triangular, or A ~ QR with Ra > 0. 

Theorem 3. Given any dQ with Q^dQ anti- symmetric, and any dR, we can compute dA = QdR-\- dQR = 
Q{dRR-^ + Q^dQ)i? = {R'^(g> Q) ((i?^®uppcr -f {Q^dQ)) . The linear maps from (lower Q'^dQ ), dR 

to dA then has determinant (nr=i CJdl^i^*) = nr=i Cr*- 



Note: It is straightforward to imagine computing the Jacobian numerically. We can make "^"^ ' perturba- 
tions to R, say by systematically adding .001 to each entry, and then compute dA via dA — QR—Q{R + dR). 

Similarly we can make n(n—l)/2 perturbation to Q by choosing two distinct columns of Q and multiplying 
the n X 2 matrix formed from these columns by a matrix such as 

cos(.GOl) sin(.OOl) " 
-sin(.GGl) cos(.OOl) ' 

The dQ is the new Q minus the old one. 

The resulting dA's may be squashed into an rt^ by matrix whose determinant could be computed 
numerically, confirming Theorem 3. 

Note: It is easy to compute the inverse map, i.e., given dA compute dQ and dR. Numerically we can ask 
for [Q, R] = qr{A) and [Q + dQ,R + dR] = qr{A + dA). 

Mathematically, dA = Q{dRR~^ + Q^dQ)R so if lowcr(M) denotes the strictly lower part of M and 
upper(M) denotes the upper part of M, then 

dQ = Q(lowcr(Q^dyli?-i) -lower(Q^dAi?-i)^) 
and dR = (upper(Q^dAi?-i) + lower(Q^dAdi?~i)^) i? 



6 Jacobians for Spherical Coordinates 

The Jacobian matrix for the transformation between spherical coordinates and Cartesian coordinates has a 
sufficiently interesting structure that we include this case as our final example. The structure is known as a 
"lower Hessenberg" structure. 

We recall that in M", we define spherical coordinates r, 6*1, 6*2, ... , On-i by 

xi — r cos 6i 

X2 ~ r sin 6i cos 62 

X3 — r sin 9i sin 62 cos 03 



^n— 1 
Xr) 



r Sinful sin6'2 ■ • • sm6'„_2 cos^n-i 
r sin 61 sin 62 ■ ■ ■ sin 0„_2 sin 



OT ioi i — 1, . . . , n, Xj may be written as 



J|sin( 

.i=l 



COS ( 



9n = 0) . 



We schematically place an x in the matrix below if Xj depends on the variable heading the columns. The 
dependency structure is then 



Xi 
X2 



r 

( X 

X 



7n-2 Pn-l 



Xn^2 




X 


X 


X 


X 




'^n— 1 




X 


X 


X 


X 


X 




[ 


X 


X 


X 


X 


X / 



Therefore the nonzero structure of the Jacobian matrix is represented by the pattern above. 

Matrices that are nonzero only in the lower triangular part and on the superdiagonal are referred to as 
lower Hessenberg matrices. 



It is fairly messy to write down the exact Jacobian matrix. Fortunately, it is unnecessary to do so. We 
obtain the LU factorization of the matrix by defining auxiliary variables yi as follows: 



yi = 

ys = 



sin^ 6i 

sin^ Qx sin^ 6*2 



sin^ 6*1 sin^ 6*2 •• • sin^ Qn-\ 



Differentiating and expressing the relationship between Ayi + dxj for i, j = 1, . . . , n in matrix form wo obtain 
that the Jacobian matrix Jx^y from x to y is triangular: 



V 



dxo 



\dXnJ 



fdy,\ 

dy2 



\dyn/ 



We recognize that we have an upper triangular matrix in front of the vector of Cartesian differentials and a 
lower triangular matrix in front of the vector of spherical coordinate differentials. 



/I 1 



1\ 



V 



1 ••• 1 



1/ 



/xi 



V 



\ 



X2 



fdxA 



dX2 



XnJ \dXnJ 



fdy,\ 
dy2 



\dynj 



Similarly the Jacobian matrix Jsphcricai^y from spherical coordinates to y is triangular 
/2r 



2r sin 9iXi 



2r sin 01 sin 62X2 



\ / dr \ /dyi\ 
dd2 



2rsin9i ■ ■ ■ sm0n-iXn-J V^^n-iJ 



dy2 

\dyn/ 



Therefore 



J-^ J 



The LU (really {L U) ) allows us to easily obtain the determinant as the ratio of the triangular 
determinants. The Jacobian determinant is then 



dix 



1 1 ■ ■ ■ T J 



d{r 



2"r"(sin6'i)"-Hsin6l2)"-^(sin6'„_i)xi ■ ■ ■ Xn~-i 

2^ X\ • • • Xfi 
r"-i(sin0i)"-2...(sin0„_2). 



In the next section we will introduce wedge products and show that they are a convenient tool for handling 
curvilinear matrix coordinates. 



7 Wedge Products 



There is a wedge product notation that can facihtate the computation of matrix Jacobians. In point of 
fact, it is never needed at all. The reader who can follow the derivation of the Jacobian in Section 5 is 
well equipped to never use wedge products. The notation also expresses the concept of volume on curved 
surfaces. 

For advanced readers who truly wish to understand exterior products from a full mathematical viewpoint, 
this chapter contains the important details. We took some pains to write a readable account of wedge 
products at the cost of straying way beyond the needs of random matrix theory. 

The steps to understanding how to use wedge products for practical calculations are very straightforward. 

1. Learn how to wedge two quantities together. This is usually understood instantly. 

2. Recognize that the wedge product is a formalism for computing determinants. This is understood 
instantly as well, and at this point many readers wonder what else is needed. 

3. Practice wedging the order n? or mn or some such entries of a matrix together. This requires mastering 
three good examples. 

4. Learn the mathematical interpretation of the wedge product of such quantities as the lower triangular 
part of Q^dQ and how this relates to the Jacobian for QR or the symmetric eigendecomposition. We 
provide a thorough explanation. 

8 The Mechanics of Wedging 

The algebra of wedge products is very easy to master. We will wedge together differentials as illustrated in 
this example: 

(2da:; + x'^dy + 5dw + 2dz) A {ydx - xdy) = 

{-2x - x'^y)dx Ady + 5y{dw A dx) - 5x{dw A dy) - 2y{dx A dz) + 2x{dy A dz) 

Formally the wedge product acts like multiplication except that it follows the anticommutative law 

(du A dv) = {~dv A du) 

Generally 

{pdu + qdv) A {rdu + sdv) — {ps — qr){du A dv) 

It therefore follows that 

dw A du = . 



In general 

fi{x) dxi A gj{x) dxj = ^^(•^^(^^fj^^) ~ fji^)M^)) da;, A dxj = ^ 



i<j i<j 

We can wedge together more than two differentials. For example 

2d2: A {3dx + 5dy) A 7(da; + dy + dz) = 70dx A dy A dz 

Let us write down concretely the wedge product. 

Rewriting the two examples above in matrix notation, we have 



fi{x) gi{x) 
fj{x) gj{x) 



dxi A dxj 



/ 2 y\ 

2 

X —X 

5 

V 2 ; 



2 3 7 

and F = I 5 7 
7 



In the first case, we have computed ah 2 x 2 subdeterminants and in the second case we compute the one 
3x3 determinant. In general, if F £ we compute all (^) subdeterminants of size p. 

If F{x) £ M"^^ and dx is the vector (dxi,da;2, . . . jdzn)"^ , then we can wedge together the elements of 
F{x)^dx. The result is 

}\{F{xrdx),= ^(i ? Ma^^.A-.-Adx,^, 

i=l ii<i2<---<ip 

where ^ (^^-^ ^2 ' ' ' ^ denotes the subdetcrminant of F obtained by taking rows ii,i2, ■ ■ ■ ,ip and 

all p columns. In almost MATLAB notation, ^ (^^^ ^2 ' ' ' ^ ^ dct(_F[ii i2 ... ip, :]). In simple 

English, if wc wedge together p differentials which we can store in the columns of F, then the wedge product 
computes all p x p subdeterminants. 
We use the notation 

[Fixfdxr = ^(Fixfdx),, 
1=1 

i.e., "( )^" denotes wedge together all the elements of the vector inside the parentheses. 

A notational note: We are inventing the "( )^" notation. Books such as [4] use only parentheses 

to indicate the wedge product of the components inside. Unfortunately, we have found that parenthesis 
notation is ambiguous especially when we want to wedge together the elements of a complicated expression. 
Once we decided on placing the wedge symbol "A" somewhere, we experimented with upper/lower left and 
upper/lower right, finally settling on the upper right notation as above. 

We will extend the "( notation from vectors to matrices of differentials. Wc will only wedge 

together the independent elements. For example 



(dM)^ = /\ dM„ M £ 



l<i<r 
l<j<n 



Wc use subscripts to indicate which elements to wedge over. For example (dS')f>j = /\i<:jdSij. When 
unnecessary as in the cases below, we omit the subscripts. 



(d^)^ = 


A ds,, 

l<i<j<n 


Sen 




symmetric 


(dA)'^ = 


A dA,, 

l<i<j<n 


Ael 


^nn 


antisymmetric 


(dA)'^ = 


A dAu 

l<i<n 


A e 1 




diagonal 


(dC/)^ = 


A dU,, 

l<i<i<n 


U e ] 


^nn 


upper triangular 


(dC/)^ = 


A dU,, 

l<i<j<n 


U e ] 


^nn 


strictly upper triangular 


etc. 











Definitional note: We have not specified the order. Therefore the wedge product of elements of a matrix is 
defined up to "+" or " sign. 

Notational note: We write (dMi)^(dAf2)^ for (dAfi)^ A (dA/2)^. 



9 Jacobians with wedge products 

9.1 Wedge Notation not useful (Y = X"^) 
Consider the 2x2 case of y = as in Example 1 of Section 2.2. 



With X = 



p q 

r s 



dY = XdX + AXX 



p q 




dp dq 




dp dq 




p q 


r s 




dr ds 


+ 


dr ds 




r s 



Thus 

dYii = 2pdp + qdr + rdq 

dY2i = rdp + {p + s)dr + rds 

dYi2 = qdp + {p + s)dq + qds 

dY22 = qdr + rdq + 2sds . 

If we wedge together all of the elements of Y, we obtain the determinant of the Jacobian: 

(dr)^ = A{p + sf{sp- qr) . 
The reader should compare the notation with that from the previous chapter: 



J 



We conclude that for these examples the notation is of little value. 



dp 


dr 


dq 


ds 




2p 


q 


r 





dYn 


r 


p + s 





r 


dY2i 


q 





p + s 


q 


dY,2 





q 


r 


2s 


dY22 



(7) 



9.2 Square Matrix = Antisymmetric + Upper Triangular 

Given any matrix M e M"-", let 

A = lower(Af) - lower(A'/)^, 
where lower(M) is the strictly lower triangular part of A/. Further let 

R = uppcr(Af ) + lower(A/)^, 

where uppcr(M) includes the diagonal. We have the decomposition 



Therefore 



M -- 




A 


+ \r 




/ drii 


dri2 - 






da2i 


dr 


dM = 




dasi 


da 




V 


da„i 


da 



(antisymmetric) + (upper triangular) 



32 



dri3 - daai 
dr23 - da32 
dr33 



da 



n3 



drin - da„i \ 
dr2„ - da„2 
dr3„ - da„3 



(8) 



We can easily see that the — da.y play no role in the wedge product: 

(dAf)^ = (lower (dA) + di?)^ = (dA)^(di?)^ . 

We feel this notation is already somewhat more mechanical than what was needed in the last chapter. 
Without wedges, we would have concentrated on the n{n + l)/2 parameters in R and the n{n — l)/2 in 
lower(A) and would have written a block two by two matrix as in Section 5.2. 

The reader should think carefully about the following sentence: The realization that the upper da^i's 
play no role in (dAI)^, is equivalent to the realization that the matrix X12 plays no role in the determinant 
of the 2x2 block matrix 

y _ X\2 

X21 



9.3 Triangular < — > Symmetric 

Let S ^ U + = ^+\\ . If we consider the map from U to upper(S'). it is easy to see that the Jacobian 
determinant is 2". 

(dS')'" = (diag(dS'))'' A (strictly-upper(dS'))'' 

= 2"(diag(dC/))'^ A (strictly-upper(dC/))'' 
= 2"(dC/)^ . 



9.4 General matrix functions 

Let Y = Y(X) be a 2 X 2 matrix, where X = \ . We have that 

Va:2i X22 ' 



dY = 



l^dxu + |^dxi2 + |^dx21 + ^dx22 l^dxn + l^dxi^ + l^dX21 + ^dx22 



dx\2 ^'^ dx2i dx22 dxii dxi2 ^'^ dx2i '^^ dx 

The wedge product of ah the elements is 

(dY)^ = dyii A dy2i A dyr2 A dj/22 = det f |^ ) daiii A da;2i A d^ia A da;22 = (det J){dXy 

\OXkl J 

In general, if y = Y{X) is an n x n matrix, then 

(dy)^ = (det J)(dX)^, 

where det J denotes the Jacobian determinant ■ 



9.5 Eigenproblem Jacobians 
Case I: S Symmetric 

If 5* e E"'" is symmetric it may be written 

S = QAQ^, Q^Q = /, A = diag(Ai, • • • , A„) . 

The columns of Q are the eigenvectors, while the Xi are the eigenvalues. 
Differentiating, 

dS = dQAQ'^+QdAQ'^+QAdQ^, 
Q'^dSQ = (Q^dQ)A - A(Q^dQ) + dA . 

Notice that Q^dSQ is a symmetric matrix of differentials with dA^ on the diagonal and qfdqj{Xi — Xj) as 
the ijjth entry in the upper triangular part. 

Therefore (Q'^dS'Q)'^ = n»<j " Aj|(dA)^(Q^dQ)^. We have from Example 3 of Section 4 that 
(Q'^dS'Q)^ = (dS-)^. In summary 

(dS)'' ^Y[\K~ X,\idA)'^iQ'^dQ)\ 

Many readers will wonder about the meaning of (Q^dQ)^. We did. We will not discuss this in length here 
but encourage you to ask us if you are curious. 



10 Handbook of Jacobians 



Some Jacobians that come up frequently are listed in this and the following sections. 



Function 


Sizes 


Jacobian 


V" R V id ^ 

I 1 > y\ l\ 

= {A®B)X 
"Kronecker Product" 


X TH X n 
A n X n 
B m X m 


{AY) = (detyl)"(detB)™ 


Y = X-^ 
"inverse" 


X n X n 


{AY) = {AetX)-^"{AX) 


Y = X'- 


X n X n 


{AY) = U,^^\\ + X,\{AX) 


Y ^X^ 


X n X n 


(dy) = n,>, Eto A^A^ (dX) 


Y = liAXB + B'^'XA) 
"Symmetric Kronecker" 


X,A,B nx n 
A, B sym 
B = A^ 


{AY)=U^<,\^^M,+X,Ah\{AX) 
{AY) = (detA)"+^(dX) 



11 Real Case; General Matrices 

In factorizations that involve a square m x m orthogonal matrix Q, and start from a thin, tall m x n matrix 
A (to > n always), H^dQ stands for the wedge product 

n n 

H^dQ=f\ /\ qfdq,. 

i=l 

Here is a stand-in for the first n columns of Q. If the differential matrix and the matrix whose transpose 
we multiply by are one and the same, the wedging is done over all of the columns. For example, 

n n 

V^dV^/\ /\ vfdv,, 



for y an n X n orthogonal matrix. 

The Jacobian of the last factorization (non-symmetric eigenvalue) is computed only for the positive- 
measure set of matrices with real eigenvalues and with a full set of eigenvectors. 



Factorization 


Matrix sizes & properties 


Parameter count 


Jacobian 


A = LU 


A, L,U n X n 


n2 = 






(dA) = 


lu 


L, U"^ lower triangular 
hi = 1, Vi = 1, n 


n{n—l) 1 n(n + l) 
2 ' 2 




n 

i = l 






A,R m X n 










A = QR 


Q m X m 


mn = 






(dA) = 


qr 


Q^Q = /„ 
R upper triangular 


{m.n "("+!)) + "("^H 


-1) 


U 

i=l 


'^'\dR)(H'^dQ} 




A,T, m X n 










A = UEV'^ 


U m X m, V n X n 


mn = 






(dA) = 


svd 




{mn 2 ) + "1 ^ 


n-l) 
2 


n n 

i<j i=i 


aY''" {dS}{H'^ dU){V'^ dV) 




V^v = /„ 










A,Q m X n 










A^QS 


S n X n 


mn = 






(dA) = 


polar 


Q^Q = /„ 
S positive definite 


{mn "("+!)) + "("^H 




i<j 


+ a,){dS}{Q'^dQ) 


A = XAX-' 


A,X,k n X n 


n2 = 






(dA) = 


nonsymm cig 


A diagonal 


(n.2 — n.) -I- n 




n(A. - 


Xjf(dA)(X^'^dX) 



12 Real Matrices; Symmetric "+" Cases 



Here we gather factorizations that are symmetric or, in the case when this is required, positive definite. AU 
matrices are square m x m. 



Factorization 


"+" 


Parameter count 


Matrix properties 


Jacobian 


A — LL' 
Choleski 


Positive definite 


n{n + l) _ 
2 

n(n + l) 
2 


L lower triangular 


(dA) = 

n 

2" n it'-'idL) 

i=l 


A = LDL' 
Idl 




n(n+l) _ 
2 

"("-1) ^ n 


L lower triangular 
la = I, Vi = 1, n 
D diagonal 


{dA) = 

n 

udr'{dL){dD) 


A = QAQ^ 
eig 




7i(n+l) 

2 ^ 
n(n—l) 1 

2 


Q^Q = /„ 
A diagonal 


{dA) = 
n |A, -A,|(dA)(Q^dQ) 



13 Real Matrices; Orthogonal Case 



This is the case of the CS decomposition. Note that the matrices Ui and Vi share parameters; this 

is due to the fact that the product of the first p columns of Ui and the first p rows of is invariant and 
determined. This is equivalent to saying that we introduce an equivalence relation of the set of pairs (Ui, Vi) 
of orthogonal matrices of size fc, by having 



Q 
li 



.Vi 



Q 
/, 



for any p x p orthogonal matrix Q. 

Since it is rather hard to integrate over a manifold of pairs of orthogonal matrices with the equivalence 
relationship mentioned above, we will use a different way of thinking about this splitting of parameters. We 



can assume that Ui contains 



fc(fc-i) 



parameters (a full set) while Vi only contains 



fc(fc-i) p(p-i) 



(having 



the first p columns already determined from the relationship with Ui). Alternatively, we could assign a full 
set of parameters to Vi and think of Ui as having the first r columns predetermined. In the following we 
choose to make the former assumption, and think of the undetermined part H of as a point in the Sticfel 
manifold Vj^k- 



Factorization 


Matrix sizes & properties 


Parameter count 


Jacobian 


n = k + j, p = k- j 
CS decomposition 


Q n X n 
Ui,Vi kxk 

U2,V2 j X j 

Ui,Vi,U2,V2 ort hogonal 

C = diag(cos(6i)), i = 1 . . . n 
S = diag(sin(6i)), i = 1 . . . n 


UiV^ has 

k{k 1) "^"-^^ 

U2,V2 have 
''•'-'^^ each 

C and S share j 


idA) = 

n siniOj - ei)sin{ei + ej)(de) 
i<j 

(Uf dUi)iU^ dU2}{H^ dH){V^ dV2} 



14 Real Matrices; Symmetric Tridiagonal 

The important factorization of this section is the eigenvalue decomposition of the symmetric tridiagonal. The 
tridiagonal is defined by its 2n — 1 parameters, n on the diagonal and n — I on the upper diagonal: 



T 



( a-n bn-1 

bn-1 fln-l bn^2 



02 bi 
bi ai J 



The matrix of eigenvectors Q can be completely determined from knowledge of its first row q~ (gi , . . . , g„) 
(n — 1 parameters) and the eigenvalues on the diagonal of A (another n parameters). Here dq represents 
integration over the (l,n) Stiefel manifold. 



Factorization 


Matrix sizes 


Parameter count 


Matrix properties 


Jacobian 


T = QKQT 
eig 


T,Q,A nx n 


2n-l = 
{n-l)+n 


A diagonal 


[dT) = {dA){dq) 

n 1^ 

i = l 



Proof by MATLAB 

At a basic level, as we discussed earlier, the Jacobian is the magnification in the volume element that occurs 
when a matrix is perturbed in "every conceivable direction". Hence the theoretical results can be verfied 
using actual MATLAB experiments. Consider the symmetric eigenvalue problem on the third row of the table 
in Section 12. The function jsymeig.m allows you to enter a real symmetric matrix A and computes the 
theoretical Jacobian and the numerical answer. You should easily be able to verify that both these answers 
match. We encourage you to verify the theoretical answers for other factorizations. 

Once we have the Jacobians computing the eigenvalue densities is fairly straightforward. The densities 
for the Hermite, Laguerre and Jacobi ensembles is tabulated here as well. 



function j = jsymeig(a) 

ZJSYMEIG Jacobian for the symmetric eigenvalue problem 

7„ J = JSYMEIG(A) returns the Numerical and Theoretical Jacobian for 

7„ the symmmetric eigenvalue problem. The numerical answer is 

7„ computed by perturbing the symmetric matrix A in the N*(N+l)/2 

7„ independent directions and computing the change in the 

% eigenvalues and eigenvectos. 

% 

% A is a REAL N x N symmetric matrix 

% Example: Valid Inputs for A are 
% 1) G = randn(5); A = (G + G') 

% 2) G = randn(5,10); A = G * G' 

7. 

7o J is a 1 X 2 row vector 

7o J(l) is the theoretical Jacobicin computed using the formula in the 
7o third row of the table in Section 12 of Handout # 6 

7. 
% 

7. References: 

7. [1] Alan Edelman, Handout 6: Essentials of Finite Random Matrix 
7. Theory, Fall 2004, Course Notes 18.338. 

7. [2] Alan Edelman, Random Matrix Eigenvalues. 

7. 

7o Alan Edelman aoid Raj Rao, Sept. 2004. 

7. iRevision: 1.1 $ $Date: 2004/09/28 17:11:18 $ 



format long 

[q,e] = eig(a) ; 7o Compute eigenvalues and eigenvectors 
e = diag(e); 

epsilon = le-7; 7o Size of perturbation 
n = length (a) ; °L Size of Matrix 
jacmatrix = zeros (n* (n+1) /2) ; 

k = 0; mask = triu( ones(n),l); mask = logical (mask( :)) ; 
for i = l:n, 

for j = i:n 

k = k+1; 

E = zeros (n) ; 

E(i,j) = 1; E(j,i) = 1; 

aa = a + epsilon * E; 

[qq,ee] = eig(aa) ; 

de= (diag(ee) -e) /epsilon; 

qdq = q'*(qq-q) /epsilon; qdq = qdq(:); 

jacmatrixCl :n,k) = de; 

jacmatrixC (n+1) : end,k) = qdq (mask) ; 

end 

end 



% Numerical answer 
j = l/det(jacmatrix) ; 



7. Theoretical answer 
[el,e2] = meshgrid(e , e) ; 
z = abs(el-e2); 

j = abs([j prod( z(mask) ) ]); 



Code 2 



Random Matrix Ensembles 

Joint Eigenvalue Distributions: c I A I fl V( Xj ) 



Technical name 


Traditional name 


P 


Field 


Property \ \ 


XXClIllllC CllOClllUlCiS 


vrctUsMctll cllocIllUlcs 








va^-e"^''' 

V ^^A- / — C 


GOE 
GUE 


1 

2 


R 

C 


Symmetric \ 
Hermitian / 




GSE 


4 


H 


Self-Dual / 


Laguerre ensembles 


Wishart ensembles 

Wishart real 


1 


R 


Positive 


{ 


a =-2-(n - m+1) - 1 


Wishart complex 
Wishart quaternion 


2 
4 


C 
H 


Semi-Definite^ 




Jacobi ensembles 


MANOVA ensembles 








V(X)= X^(l-X)'' 


MANOVA real 


1 


R 


Positive 




a =-|-(nj - m+1) - 1 


MANOVA complex 


2 


C 


Semi-Definite 




b =-^(^2 ~ i^+l) ~ 1 


MANOVA quaternion 


4 


H 





Invariance 



QTaQ 

u"au 

S^AS 





A^ QTaQ 
A^ U"AU 
A^S°AS 



Connected Matrix Problem 



(EIG) 

Symmetric 
Eigenvalue Problem 



(SVD) 

Singular Value 
Decomposition 



(QZ) 

Generalized Symmetric 
Eigenvalue Problem 



Type of ensemble 



MATLAB code 



Hermite 



A = randn(n); A = (A+A')/2; 

A = randn(n) + i * randn(n); A = (A+A') / 2; 

X = randn(n) + i * randn(n); Y = randn(n) + i * randn(n); 

A=[X Y; -conj(Y) conj(X)]; A = (A+A') 7 2; 



Laguerre 



A = randn(m, n); A = A*A'; 
A = randn(m, n) + i * randn(m, n); 
X = randn(m, n) + i * randn(m, n); 
A=[X Y; -conj(Y) conj(X)]; 



A = A* A'; 

Y = randn(m, n) + i * randn(m, n); 
A = A* A'; 



Jacobi 



X = randn(m, n j); Y = randn(m, n j); A = (X * X') / (X * X' + Y * Y'); 



Y = randn(m, n, ) + i * randn(m, n,); A = (X * X') / (X * X' + Y * Y'); 



X = randn(m, n j) + i * randn(m, n j ) 
Xj = randn(m, n j) + i * randn(m, n ^ ) 
Yj = randn(m, n^) + i * randn(m, n^) 
X=[X, X^; -conjCX^) conj(Xj)]; Y = [Y, Y^; -conj(Y^) conj(Y,)]; 
A = (X * X') / (X * X' + Y * Y'); 



X 2= randn(m, ) + i * randn(m, n j); 
Y 2= randn(m, ) + i * randn(m, n 



15 Volume and Integration 



15.1 Integration Using Differential Forms 

One nice property of our differential form notation is that if y ~ y{x) is some function from (a subset of) 
to M", then the formula for changing the volume element is built into the identity 

/ /(y)dj/i A... Ad?/„ = / f{y{x))dxl^...^Axn, 

Jy(S) JS 

because as we saw in 9.4 the Jacobian emerges when we write the exterior product of the dj/'s in terms of 
the Ax's. 

We will only concern ourselves with integration of n- forms on manifolds of dimension n. In fact, most 
of our manifolds will be flat (subsets of M"), or surfaces only slightly more complicated than spheres. For 
example, the Stiefel manifold Kn,n of n by p orthogonal matrices Q {Q^Q ~ Im) which we shall introduce 
shortly. Exterior products will give us the correct volume element for integration. 

If the Xi are Cartesian coordinates in n-dimensional Euclidean space, then (dx) = dxi A dx2 A . . . dx„ is 
the correct volume element. For simplicity, this may be written as da;ida;2 . . . da;„ so as to correspond to the 
Lebesgue measure. Let qi be the ith component of a unit vector q G M". Evidently, n parameters is one too 
many for specifying points on the sphere. Unless g„ = 0, we may use qi through qn-i as local coordinates 
on the sphere, and then dqn may be thought of as a linear combination of the dqi for i < n. qidqi = 
because q'^q = 1). However, the Cartesian volume element dqidq2 . . .dqn-i is not correct for integrating 
functions on the sphere. It is as if we took a map of the Earth and used Latitude and Longitude as Cartesian 
coordinates, and then tried to make some presumption about the size of Greenland^. 

Integration: 

Ixes f(^)i'^-'^) '3'' Is fi'^-^) other related expressions will denote the "ordinary" integral over a region 

5" el. 

Example. exp(-||a;||V2)(dx) = (27r)"/2 and similarly /g„,„ exp(-||a;|||,/2)(dyl) = (27r)"'/2. = 
ti{A'^A) = J2i j '^ij ~ "Frobenius norm" of A squared. 

If an object has n parameters, the correct differential form for the volume element is an 7i-form. What 
about X e S^'-S i.e., {x e M" : ||.t|| = 1}? Ai=i dx^ = (dx)'^ = 0. We have J^^i = 1 ^ E xtdxi = ^ 
dxn — —-^{xidxi + ■ ■ ■ + Xn-idxn^i). Whatever the correct volume element for a sphere is, it is not (dx). 

As an example, we revisit spherical coordinates in the next section. 

15.2 Plucker Coordinates and Volume Measurement 

Let F £ W-'P. We might think of the columns of F as the edges of a parallelopipcd. By defining P1(F) 
("Plucker(F)" ), we can obtain simple formulas for volumes. 

Definition 1. Pl(-F) is the vector of p x p subdeterminants of F written in natural order. 

do not think that I have ever seen a map of the Earth that uses Latitude and Longitude as Cartesian coordinates. The 
most familiar map, the Mercator map, takes a stereographic projection of the Earth onto the (complex) plane, and then takes 
the image of the entire plane into an infinite strip by taking the complex logarithm. 
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Definition 2. Lei vol(_F) denote the volume of the parallelopiped {Fx : < Xi < 1}, i.e., f/ie volume of the 
parallelopiped with edges equal to the columns of F. 

Theorem 4. vol(F) = IliLi "'i ~ det(F-^F)-'^/^ = IliLi '^j* ~ I|P1(^)II ' where the Ui are the singular values 
of F , and the ru are the diagonal elements of R in F ~ Y R, where Y G R"'^ has orthonormal columns and 
R is upper triangular. 

Proof. Let F = UY.V^ , where U G K"^^ and V G IRp-p has orthonormal columns. The matrix a denotes 
a box with edge sides ai so has volume Jlcr;. The volume is invariant under rotations. The other formulas 
follow trivially except perhaps the last which follows from the Cauchy-Binet theorem taking A = and 
B = F. □ 

Theorem 5. (Cauchy-Binet) Let C = AB be a matrix product of any kind. Let M(^J'*'' ) denote the px p 
minor 

det{Mi^j^)i<k<p,i<i<p- 

Ln other words, it is the determinant of the submatrix of M formed from rows ii, . . . ,ip and columns ji, . . . ,jp. 



The Cauchy-Binet Theorem states that 



^ I *i J ■ ■ ■ J *p 
fci , . . . , /Cp 



E 



A 



3l<3-2<---<3p 



. . . ^Ip 
Jl I ■ • ■ I Jp 



B 



111 
fci, 



7 Jp 
k 



Notice that when p = \ this is the familiar formula for matrix multiplication. When all matrices are p x p, 
then the formula states that 

det C = det A det B . 

Corollary 6. Let F e IR"^^ have orthonormal columns, i.e., F^F = Ip. Let X £ IR"-^'. // span(F) = 
span(X), then vol(X) = det(F^X) = P1(F)^ • P1(X). 

Theorem 7. Let F,V e M"'^ be arbitrary. Then 

p 

P1(J^)'^P1(F) = det(F'^X) = |vol(i^)| |vol(F)|J|cos6'i, 

i=l 

where 61, . . . ,9p are the principal angles between span{F) and span(V). 

Remark 1. Lf Sp is some p- dimensional surface it is convenient for F^'^'^ to be a set of p orthonormal tangent 
vectors on the surface at some point a;*^*^ and V^'^ to be any "little" parallelopiped on the surface. 



If we decompose the surface into parallelopipeds we have 

vol(5p) « ^Pl(i^('))^Pl(yW) . 

and 

J /(x)d(surface) « ^ Pl(F«)^Pl(y(')) = ^ dct((i^«)^y«) 

Mathematicians write the continuous hmit of the above equation as 

J /(x)d(surface) = ^ /(j;)(F^da;)^ . 
Notice that (F(x)"^da;)'^ formally computes Pl(i^(a::)). Indeed 



dxi^ A ... A dxip 



ii<...<ip 



III. Generalized dot products algebraically: Linear functions l{x) for x G M" may be written l{x) = J2 fi^i^ 
i.e., f'^x. 

The information that specifies the linear function is exactly the components of /. 

Similarly, consider functions 1{V) for V € R"'^ that have the form 1{V) = dct(i^^), where F e 
The reader should check how much information is needed to "know" 1{V) uniquely. (It is "less than" 
all the elements of F.) In fact, if we define 

Eij^i2...ip = The n X p matrix consisting of columns ii, . . . , of the identity, 

i.e., {E ~ eye(7i); i?i^i2^...,ip = E{:,[ii . . .ip])) in pseudo-Matlab notation, then knowing l{Ei-^^,,,^i^) 
for all zi < ... < is equivalent to knowing I everywhere. This information is precisely all p x p 
subdeterminants of F, the Plucker coordinates. 

IV. Generalized dot product geometrically. If F and V E E"'P we can generalize the familiar f'^v = 
11/11 \\v\\cos9. 

It becomes 

p 

Pl(i^)^Pl(y) = det(F'^y) = |vol(F)| |vol(y)| J|cos6'i . 

2=1 

Here vol(A/) denotes the volume of the parallclopiped that is the image of the unit cube under M (Box 
with edges equal to columns of A/). Everyone knows that there is an angle between any two vectors 
in R". Less well known is that there arc p principal angles between two p-dimensional hypcrplanes in 
R". The 9i arc principal angles between the two hypcrplanes spanned by F and by V. 

Easy special case: Suppose F'^F = Ip and span(F) = span(y). In other words, the columns of F are 
an orthonormal basis for the columns of V. In this case 

det{F^V) = vol(V) . 

Other formulas for yol(V) which we mention for completeness is 

p 

i=l 

the product of the singular values of V. And 

vol(y) 



\ Ji<...- 



1...P 

ii . . .ip 



■ <ip 

the sum of the p x p subdeterminants squared. In other words vol(V^) = ||Plucker(y)|| 



1. Volume measuring function: (Volume element of Integration) 

We now more explicitly consider V g as a parallelopiped generated by its columns. For general 
F G det{F'h') may be thought of as a "skewed" volume. 

We now carefully interpret /\^^i{F^<ix)i, where F £ and dx = (dxida::2 . . ■dxn)'^ ■ We remind the 
reader that {F'^dx)i = X]j'=i Fji^Xj. 

By /\!l^i{F^dx)i we are thinking of the linear map that takes V g E"'P to the scalar det(i^"^). In our 
mind we imagine the following assumptions 

• We have some p-dimensional surface Sp in R" . 

• At some point P on Sp, the columns of X e M"'P are tangent to iSj,. We think of X as generating a 
little parallelopiped on the surface of Sp. 

• We imagine the columns of F are an orthonormal basis for the tangents to Sp at P. 

• Then /\i^i{F^dx)i — (F^dx)^ is a function that replaces parallclopipcds on Sp with its Euclidean 
volume. 

Now wc can do "the calculus thing." Suppose we want to compute 




/(x)d(surface) , 



the integral of some scalar function / on the surface Sp. 

We discretize by circumscribing the surface with little parallclopipcds that we specify by using the matrix 
y(0 g K»,P. The word circumscribing indicates that the p columns of V^"^"^ are tangent (or reasonably close 
to tangent) at a:^*'' . 

Let be an orthonormal basis for the tangents at a;*^*-* . Then / is discretized by 

/w^/(x('))det(i^(')VW) 

i 

We then use the notation 

f I{x)l\{F^dx\= I f{x){F^dxr 

to indicate this continuous limit. Here / is a function of x. 

The notation does not require that F be orthonormal or tangent vectors. These conditions guarantee 
that you get the correct Euclidean volume. Let them go and you obtain some warped or weighted volume. 
Linear combinations are allowed too. 

The integral notation does require that you feed in tangent vectors if you discretize. Careful mathematics 
shows that for the cases we care about, no matter how one discretizes, the limit of small parallclopipcds 
gives a unique answer. (Analog of no matter how you take small rectangles to compute Ricmann integrals, 
in the limit there is only one unique area under a curve.) 

15.3 Overview of special surfaces 

Wc arc very interested in the following three mathematical objects: 

1. The sphere = {x : \\x\\ = 1} in M" 

2. The orthogonal group 0{n) of orthogonal matrices Q {Q^Q = I) in M"". 

3. The Sticfcl manifold of tall skinny matrices Y G IR"^ with orthogonal columns (y-^ = Ip). 



Of course the sphere and the orthogonal group are special cases of the Stiefel manifold with p = 1 and p = n 
respectively. 

We will derive and explain the meaning of the following volume forms on these objects 

Sphere {H'^dq)^^2 n where H is a. so-called Householder transform 

Orthogonal Group (Q^dQ)^ = (Q^dO)f>j. 

Stiefel Manifold (Q^dy)f>j where Ip. 



The concept of volume or surface area is so familiar that the reader might not see the need for a formalism 
to encode this concept. 

Here are some examples. 



Xi 
X2 



L 



2^11 = 1 



Example 1: Integration about a circle Let x = 

f{x){~xidxi + X2dx2) = / f{cos9,sm0)d6 

Jo 

Algebraic proof: ) ( sin^ ) ^P^^*^*^ 

/ ^X2 \ / dx] \ ( — sin0 \ ^ ( — sin0 , ,^ 
-.,d.,+xrd.,^d^= :\[^:\ = [ eosJ cos^ 



Geometric Interpretation: Approximate the circle by a polygon with k sides. 

^/(x(')) \^ (xW-a;('-i))« / /(cos0,sin0)d0. 



Note that the dot product computes 

Geometric interpretation of the wrong answer: Why not J j(x)dx\l 

\\x\\=\ 

This is approximated by 

T 

; 

which is the integral of / over the "shadow" of the circle on the x-axis. 

Example 2: Integration over a sphere Given q with IjgH = 1, let i?((7) be any n x n orthogonal matrix 
with first column q. (There are many ways to do this. One way is to construct the "Householder" reflector 

T 

H{q) = I — 2^J|rj, where v = ei — q, and ei is the first column of /.) 
The sphere is an n — 1 dimensional surface in n dimensional space. 
Integration over the sphere is then 



where dS is the surface "volume" on the sphere. 

Consider the familiar sphere n — 3. Approximate the sphere by a polyhedron of parallelograms with k 
faces. A parallelogram may be algebraically encoded by the two vectors forming edges, into a matrix V{i). 
We have ^ f{x{vi)) ■ det[H{q)'^) approximates this integral. 

Example 3: The orthogonal group Let 0{n) denote the set of orthogonal n x n matrices. This is an 

2 

dimensional set living in M" . 



Tangents consist of matrices T such that Q^T is anti-symmetric. We can take an orthogonal, but not 
orthonormal set to consist of matrices Q{Eij — Eji) for i > j {{Eij^^ = 1) ii i = k and j = I; otherwise). 
The matrix "F'V roughly amounts to taking twice the triangular part of Q^dQ. 

Let 0{m) denote the "orthogonal group" of to x m matrices Q such that Q^Q = I. We have seen 
(Q^dQ) = Ai>j iT'^lj is the natural volume element on 0{M). Also, notice that 0{n) has two connected 
components. When to = 2, we may take 




for half of 0(2). This gives 

n'^dO-f sm9\f-sm9d9 -cos6'd6'\ / -dO \ 

^ ^~y-sm9 cos9 )y cos9d9 -sinM^ ) ^ \ d9 ) 

in terms of d9. Therefore (Q^dQ)^ = d9. 



15.4 The Sphere 



Students who have seen any integral on the sphere before probably have worked with traditional spherical 
coordinates or integrated with respect to something labeled "the surface element of the sphere." We mention 
certain problems with these notations. Before we do, we mention that the sphere is so symmetric and so 
easy, that these problems never manifest themselves very seriously on the sphere, but they become more 
serious on more complicated surfaces. 

The first problem concerns spherical coordinates: the angles are not symmetric. 

They do not interchange nicely. Often one wants to preserve the spherical symmetry by writing x = qr, 
where r = \\x\\ and g = Of course, q then has n components expressing n — 1 parameters. The 

n quantities dqi, . . . ,dqn are linearly dependent. Indeed differentiating q^'^q = 1 we obtain that q^dq = 
Er=i = 0. 

Writing the Jacobian from x to q,r is slightly awkward. One choice is to write the radial and angular 
parts separately. Since da; = qdr + dgr. 



■j'^dx 



dr 



and (/ — qq'^)dx ~ rdq . 



We then have that 



da; = dr A (rdq) 



-Mr(dg), 



where (dq) is the surface element of the sphere. 

We introduce an explicit formula for the surface element of the sphere. Many readers will wonder why 
this is necessary. Experience has shown that one can need only set up a notation such as dS* for the surface 
element of the sphere, and most integrals work out just fine. We have two reason for introducing this formula, 
both pedagogical. The first is to understand wedge products on a curved space in general. The sphere is 
one of the nicest curved spaces to begin working with. Our second reason, is that when we work on spaces 
of orthogonal matrices, both square and rectangular, then it becomes more important to keep track of the 
correct volume element. The sphere is an important stepping stone for this understanding. 

We can derive an expression for the surface element on a sphere. We introduce an orthogonal and 
symmetric matrix H such that Hq = ei and Hei = q, where ei is the first column of the identity. Then 



Hdx — e\dr + Hdqr = 



/ dr \ 

riHdq)2 
r{Hdq)3 



\r{Hdq)nJ 



Thus 

n 

(da;)'" = (Fdx)'" = r"-Mr /\(iJdg), . 

i=2 

We can conclude that the surface element on the sphere is (dg) = A"=2(^'^9)i- 




Householder Reflectors 

We can explicitly construct an H as described above so as to be the Householder reflector. Notice that 
H serves as a rotating coordinate system. It is critical that H{:,2 : n) is an orthogonal basis of tangents. 
Choose w = ei — g the external angle bisector and 

T 

vv 

H = I-2^. 

See Figure 1 which illustrate that H is a. reflection through the internal angle bisector of g + ei. 

Notice that {Hdq)i = and every other component J2]=i ^ij'^lj (* 7^ 1) be thought of as a tangent 
on the sphere. H = , Hq = ei,Hei = gi, -ff^ = /, and H is orthogonal. 



plane 
of 




reflection 

Figure 1: 

I think many people do not really understand the meaning of {Q'^dQ)^. I like to build a nice cube on 
the surface of the orthogonal group first. Then we will connect the notations. First of all 0{n) is an "^"^""'"•^ 
dimensional surface in ]R" . At any point Q, tangents have the form Q ■ A, where A is anti-symmetric. 



The dot product of two "vectors" (matrices) X and Y is X - Y = tr(X^) = 'Yl-^ij^ij- Q orthogonal 
QX ■ QY = X Y. 



If Q = / then the matrices Aij — {Eij — Eji)j\f2 for i < j clearly form an 



— yj^tj — j^jij/ V ^ ' J >-i>-o,i±j. o,i± "("^^^^ dimensional cube 

tangent to 0{n) at /, i.e., they form an orthonormal basis for the tangent space. Similarly the n{n — l)/2 
matrices Q[Eij — Eji)/^/2 form such a basis at Q. 

One can form an F whose columns are vec{Q{Eij — Eji)) for i < j. The matrix would be by n{n— 1)/2. 
Then F'^{dqij)i^j would be the Euclidean volume element. Mathematicians like to throw away the so 
that the volume element is off by the factor 2"("~i)/'*. This docs not seem to bother anyone. In non-vec 
format, this is {Q^dQ)^ . 



Application 

Surface Area of Sphere Computation 

We directly use the formula (dx)'^ = r'^~'^dr{Hdq)^: 



(27r)' 



e 2 



/ r"-ie-5'^'dr / (HdqY 

Jr=0 J 



2^T g) / (Hdqr 



r(|) 



and An is the surface of the sphere of radius 1. For example, 

A2 = 2tt (circumference of a circle) 
^3 = 

= 



4tt (surface area of a sphere in 3d) 
2^2 



15.5 The Stiefel Manifold 
QR Decomposition 

Let A e K"^'"(?Ti < n). Taking x — a^, we see 3Hi such that Hi A has the form 







then construct an H2 = 



(i ... 0\ 





Ho 



( 



V 



so that H2H1A = 



fx X \ 

X 



V y 







or A={Hi--- H„ 



R 



O 



We can 



V / 

Continuing i7,„ . . . Hi A 



where i? is to x m upper triangular (with positive diagonals). 



Vo ... 0/ 

let Q = the first m columns oi Hi - ■ ■ Hm ■ Then A = QR as desired. 



The Stiefel Manifold 

The set of Q e K"^^ such that Q^Q = Ip is denoted V^_„ and is known as the Stiefel manifold. Considering 
the Householder construction, 3Hi, ■ ■ ■ , Hp such that. 



1 







HpHp^i ■ ■ ■ HiQ 



so that Q ~ 1st p columns of H1H2 ■ ■ ■ Hp^iHp. 

Corollary 8. The Stiefel manifold may be specified by (n ~ I) + (n — 2) + ■ ■ ■ + {n — p) = pn — \/2p{p + 1) 
parameters. You may think of this as the pn arbitrary parameters in annxp matrix reduced by the p{p+l)/2 
conditions that qfqj — Sij for i > j . You might also think of this as 

dim{Q} =dim{A} - dim{i?} 

T T 

pn p{p + l)/2 

It is no coincidence that it is more economical in numerical computations to store the Householder parameters 
than to compute out the Q. 

This is the prescription that we would like to follow for the QR decomposition for the n by p matrix A. 
If Q G M^'" is orthogonal, let H be an orthogonal p by p matrix such that H^Q is the first p columns of /. 
Actually H may be constructed by applying the Householder process to Q. Notice that Q is simply the first 
p columns of H. 

As we proceed with the general case, notice how this generalizes the situation when p = 1. If ^ = QR, 
then d^ = QdR + dQR and H'^dA ^ H^QdR + H^dQR. Let H = [hi, . . . , /i„]. The matrix H^QdR 
is an n by p upper triangular matrix. While H^dQ is (rectangularly) antisymmetric, {hfhj = implies 
hfdhj ^ -hjdhi) 

Haar Measure and Volume of the Stiefel Manifold 

It is evident that the volume element in mn dimensional space decouples into a term due to the upper 
triangular component and a term due to the orthogonal matrix component. The differential form 

m n 

{H^dQ)=f\ /\ hfdh, 

j = l i=j+l 

is a natural volume element on the Stiefel manifold. 
We may define 

= / {H^dQ). 
Js 

This represent the surface area (volume) of the region S on the Stiefel manifold. This "measure" ^ is known 
as Haar measure when m ~ n. It is invariant under orthogonal rotations. 

Exercise. LetTm{a) = tt™'^™^-'^''/'' r[a — ^(i — !)]• Show that the volume ofVm,n is 

'2m—.rnn/2 

Vol (K^,„) = -FTTT ■ 



Exercise. What is this expression when n = 2? Why is this number twice what you might have thought 
it would be? 

Exercise. Let A have independent elements from a standard normal distribution. Prove that Q and R 
are independent, and that Q is distributed uniformly with respect to Haar measure. How are the elements on 
the strictly upper triangular part of R distributed. How are the diagonal elements of R distributed? Interpret 
the QR algorithm in a statistical sense. (This may be explained in further detail in class). 

Readers who may never have taken a course in advanced measure theory might enjoy a loose general 
description of Haar measure. If G is any group, then we can define the map on ordered pairs that sends 
(g, /i) to g~^h. If G is also a manifold (or some kind of Hausdorff topological space), and if this map is 
continuous, we have a topological group. An additional assumption one might like is that every g G G has 



an open neighborhood whose closure is compact. This is a locally compact topological group. The set of 
square nonsingular matrices or the set of orthogonal matrices are good examples. A measure n{E) is some 
sort of volume defined on E which may be thought of as nice ( "measurable" ) subsets of G. The measure 
is a Haar measure if ^i{gE) = /i(-E), for every g Cz G. In the example of orthogonal n by n matrices, the 
condition that our measure be Haar is that 

/ HQm^dQ)"^ - / f{QoQ){Q^dQ)\ 

In other words, Haar measure is symmetrical, no matter how we rotate our sets, we get the same answer. 
The general theorem is that on every locally compact topological group, there exists a Haar measure /u. 

16 Exterior Products (The Algebra) 

Let V be an ?i-dimensional vector space over R. For p = 0, 1, . . . , n we define the pth exterior product. For 
p = it is M and for p = 1 it is For p = 2, it consists of formal sums of the form 

ai{ui A Wi), 

i 

where e M and Ui,Wi G V. (We say "ui wedge w^.") Additionally, we require that {au + w) A w = 
a{u A w) + (?J A If), u A {hv + w) = b{u Av) + {u Aw) and u A u=0. A consequence of the last relation is that 
u A w — ~w A u which we have referred to previously as the anti-commutative law. We further require that 
if ei, . . . , e„ constitute a basis for V, then A Cj for i < j, constitute a basis for the second exterior product. 
Proceeding analogously, if the form a basis for V we can produce formal sums 

^ ^ c-y(e^j A A ■ ■ ■ A e^^), 

7 

where 7 is the multi- index (71, . . . ,7p), where 71 < • • • < 7p. The expression is multilinear, and the signs 
change if we transpose any two elements. 

The table below lists the exterior products of a vector space V = {ciCi}. 



p 


pth Exterior Product 


Dimension 




1 

2 
3 


y" = M 

V^ = V = {c,ej 
= {E.<j<fe Cyfcei A e.j A Cfe} 


1 
n 

n{n - l)/2 
n{n~ l)(n- 2)/6 


P 


= {T.^,<^■,<...<i^ Ci,i2...^,ei, A 6,, A ... A 6,^ } 


6 


n 

n + 1 


y" = {cei A 62 A ... A e„} 
= {0} 


1 





In this book V ~ {^Cidx;}, i.e. the 1-forms. Then consists of the p-forms, i.e. the rank p exterior 
differential forms. 
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